1 Etude de la répétabilité spectrale d’un prétraitement

#============ PARAMETRES A MODIFIER quand vous utilisez MBXUCL
## A. Paramètres globaux liés aux noms et chemin d'accès des dataset
dataname <- "HumanSerum" # nom du jeu de données
## Mettre ici le nom du chemin d'accès du répertoire où sont les spectres prétraités
# spectres prétraités
data.path <- "/Users/manon/Documents/Conférences/RFMF2018/Atelier_PepsNMR/RMD"
# load le fichier .RData créé précédemment
load(file.path(data.path,paste0(dataname,".RData")))
dataset <- Re_Spectrum_data
# ou
# dataset <- read.table(file = file.path(data.path,paste0(dataname,"_Spectra.csv")),
# sep=";", dec=" ", row.names = 1, header = TRUE,
# check.names =FALSE)
dataset <- as.matrix(dataset)
## B. Définition des groupes de spectres par sujet.
# group_HS est loadé avec le fichier .RData créé précédemment
# ou
# dataGroup.path <- system.file("extdata", package = "PepsNMR") # données du package PepsNMR
# dir(dataGroup.path) # ce que contient le répertoire data.path
# group_HS <- read.csv(file.path(dataGroup.path, "Group_HS.csv"), header = TRUE, sep = ";")
group <- as.factor(substr(group_HS[,1],4,5)) # donneur comme groupe
table(group) # fréquences de chaque groupe
## group
## D1 D2 D3 D4 
##  8  8  8  8
group_num <- group_HS[,2] # Version numerique des noms de groupes
table(group_num) # fréquences de chaque groupe
## group_num
## 1 2 3 4 
## 8 8 8 8
ls() # éléments dans l'environnement global
## [1] "data.path"        "dataname"         "dataset"         
## [4] "Fid_info"         "group"            "group_HS"        
## [7] "group_num"        "Re_Spectrum_data"
## C. Setup des paramètres de représentation graphique des paramètres de prétraitement
num.stacked <- 2 # nombre de graphes par sortie graphique
## D. paramètres des méthodes incluses dans le MIC à modifier
ncompPCA = 4 # nombre de composantes en PCA
nClust = 4 # nombre de clusters pour le clustering
nLVPLSDA = 4 # nombre de variables latentes en PLS-DA

1.1 Inertia

# INERTIA
Inertia.res = MBXUCL::Inertia(x = dataset, y = group_num, print = FALSE)
colnames(Inertia.res[["Between_within"]]) <- c("Inertia Between groups","Inertia Within groups", "Total inertia")
pander(Inertia.res[["Between_within"]])
  Inertia Between groups Inertia Within groups Total inertia
Value 6457 711.9 7169
Percentage 90.07 9.93 100
pander(Inertia.res[["Per_group"]])
  N Inertia_group Inertia_group100 Inertia_moy_group
Group 1 8 169.8 23.85 21.22
Group 2 8 256.2 36 32.03
Group 3 8 91.47 12.85 11.43
Group 4 8 194.4 27.3 24.3
Total 32 711.9 100 NA

1.2 PCA

# PCA
PCA.res = MBXUCL::SVDforPCA(dataset, ncomp=ncompPCA)

Eigenvalues:

eigval_mat <- rbind(PCA.res$eigval, PCA.res$var, PCA.res$cumvar)
rownames(eigval_mat) <- c("eigen values", "percentage of variance", "cumulated percentage of variance")
eigval_mat <- round(eigval_mat,2)
pander(eigval_mat[,1:4], caption = "valeurs propres PCA pour les 4 premières CP")
valeurs propres PCA pour les 4 premières CP
  PC1 PC2 PC3 PC4
eigen values 5067 1389 330.4 93.69
percentage of variance 70.68 19.38 4.61 1.31
cumulated percentage of variance 70.68 90.06 94.67 95.98
DrawScores(PCA.res, drawNames=TRUE, type.obj = "PCA",
createWindow=FALSE, main = paste0("PCA score plot for ", dataname),
color = group, pch = group,axes =c(1,2))

DrawLoadings(PCA.res,type.obj = "PCA",
createWindow=FALSE, main = paste0("PCA loadings plot for ", dataname),
axes = c(1:2), loadingstype="s", num.stacked = num.stacked, xlab="ppm")
## [[1]]

1.3 Unsupervised clustering

ClustMIC.res = MBXUCL::ClustMIC(Intensities = dataset, nClust = nClust,
Trcl = group_num, Dendr = FALSE)
ClustMIC.res <- ClustMIC.res[c(1,3,5,7)] # avec clustering hiérarchique
names(ClustMIC.res) <- c("Dunn index", "Davies-Bouldin index", "Rand index", "Adjusted Rand index")
pander(ClustMIC.res, caption ="Clustering hiérarchique")
Clustering hiérarchique
Dunn index Davies-Bouldin index Rand index Adjusted Rand index
0.8262 0.6712 1 1

1.4 PLS-DA

PLSDA.res = PLSDA(x = dataset, y = group_num, nLV = nLVPLSDA, drawRMSEP = TRUE)

perf.plsda = PLSDA.res[4:6]
perf.plsda <- matrix(unlist(perf.plsda) , ncol=4, byrow = TRUE,
dimnames = list(names(perf.plsda), names(perf.plsda[[1]])))
perf.plsda <- round(perf.plsda,2)
pander(perf.plsda, caption = "Validation PLS-DA")
Validation PLS-DA
  Y1 Y2 Y3 Y4
RMSEP 0.1 0.08 0.1 0.11
R2 0.98 0.99 0.98 0.97
Q2 0.94 0.96 0.94 0.93
par(mfrow=c(1,1),xpd=TRUE)
barplot(perf.plsda[1,], beside = TRUE, col="blue", main = "RMSEP")

barplot(perf.plsda[2:3,], beside = TRUE, col=c("blue","red"), main = "R2 and Q2")
legend("bottomright", legend = c("R2", "Q2"), col = c("blue","red"),
pch = 15)

pander("Q2 cumulé pour toutes les réponses")

Q2 cumulé pour toutes les réponses

mean(perf.plsda["Q2",])
## [1] 0.9425
DrawScores(PLSDA.res, drawNames=TRUE, type.obj = "PLSDA",
createWindow=FALSE, main = paste0("PLS-DA score plot for ", dataname),
color = group, pch = group,axes =c(1,2))

DrawLoadings(PLSDA.res, type.obj = "PLSDA",
createWindow=FALSE, main = paste0("PLSDA loadings plot for ", dataname),
axes = c(1:2), loadingstype="s", num.stacked = num.stacked, xlab = "ppm")
## [[1]]

2 Evaluation de la répétabilité de différentes corrections de la Baseline

#============ PARAMETRES A MODIFIER quand vous utilisez MBXUCL
## A. Paramètres globaux liés aux noms et chemin d'accès des dataset
dataname <- "HumanSerum" # nom du jeu de données
## Mettre ici le nom du chemin d'accès du répertoire où sont les spectres prétraités
## Mettre ici le nom du chemin d'accès du répertoire où sont les FID
data.path <- system.file("extdata", package = "PepsNMR") # données du package PepsNMR
dir(data.path) # ce que contient le répertoire data.path
## [1] "Group_HS.csv" "HumanSerum"
# ## B. Définition des groupes de spectres par sujet.
# group_HS <- read.csv(file.path(data.path, "Group_HS.csv"), header = TRUE, sep = ";")
# group <- as.factor(substr(group_HS[,1],4,5)) # donneur comme groupe
# table(group) # fréquences de chaque groupe
# group_num <- group_HS[,2] # Version numerique des noms de groupes
# table(group_num) # fréquences de chaque groupe
## C. Setup des paramètres de représentation graphique des paramètres de prétraitement
num.stacked <- 2 # nombre de graphes par sortie graphique
## D. paramètres des méthodes incluses dans le MIC à modifier
ncompPCA = 4 # nombre de composantes en PCA
nClust = 4 # nombre de clusters pour le clustering
nLVPLSDA = 4 # nombre de variables latentes en PLS-DA

2.1 Fonction de prétraitement depuis la correction de la baseline

# Paramètres de la baseline qu'on veut faire varier
all_cond_BC <- expand.grid(lambda = c(1e5, 1e6,1e7,1e8), p = c(0.1,0.05, 0.01))
pander(all_cond_BC, caption = "Paramètres de correction de la baseline à faire varier")
BCconditions <- paste(paste("L:", all_cond_BC$lambda), paste("p:", all_cond_BC$p), sep=" ")
## PepsNMR prétraitement
# Lecture des Fids
fidList <- ReadFids(path = file.path(data.path, dataname), subdirs = FALSE)
# créer la matrice spectrale
Fid_data <- fidList[["Fid_data"]]
# créer la matrice d'infos sur les paramètres d'acquisition
Fid_info <- fidList[["Fid_info"]]
# réarranger les noms de ligne dans Fid_data et Fid_info pour les
# rendre plus lisibles
rownames(Fid_data) <- substr(rownames(Fid_data),1,5)
rownames(Fid_data) <- sub("-","", rownames(Fid_data))
rownames(Fid_info) <- rownames(Fid_data)
# GroupDelayCorrection
Fid_data_GDC <- GroupDelayCorrection(Fid_data, Fid_info)
# Solvent Suppression
Fid_data_SS <- SolventSuppression(Fid_data_GDC, returnSolvent = FALSE)
# Apodization
Fid_data_Apod <- Apodization(Fid_data_SS, Fid_info)
# FourierTransform
Spectrum_data_FT <- FourierTransform(Fid_data_Apod, Fid_info)
# ZeroOrderPhaseCorrection
Spectrum_data_ZOPC <- ZeroOrderPhaseCorrection(Spectrum_data_FT)
# InternalReferencing
Spectrum_data_IR <- InternalReferencing(Spectrum_data_ZOPC, Fid_info)
#============== After Alignement pre-processing
AfterAlign_preprocess <- function(lambda.bc , p.bc) {
## Baseline correction avec les paramètres lambda.bc et p.bc à faire varier
Spectrum_data <- BaselineCorrection(Spectrum_data_IR, lambda.bc = lambda.bc,
p.bc = p.bc)
## Suppression des valeurs négatives
Spectrum_data <- NegativeValuesZeroing(Spectrum_data)
## Alignement - Warping
Spectrum_data <- Warping(Spectrum_data,reference.choice = "before")
## WindowSelection
Spectrum_data <- WindowSelection(Spectrum_data, from.ws = 10, to.ws = 0.2)
# Bucketing avec le Window selection intégré
Spectrum_data <- Bucketing(Spectrum_data, intmeth = "t", mb=500)
## Suppression des régions non informatives
Spectrum_data <- RegionRemoval(Spectrum_data, typeofspectra ="serum")
## Normalisation
Spectrum_data <- Normalization(Spectrum_data, type.norm="mean")
# Exportation des résultats
Re_Spectrum_data=Re(Spectrum_data)
return(Re_Spectrum_data)
}
# application de la fonction AfterAlign_preprocess sur les spectres
# prétraités jusqu'à InternalReferencing avec les différents parmaètres de
# la correction de la baseline
res_Preprocessing <- mapply(FUN = AfterAlign_preprocess, lambda.bc=all_cond_BC$lambda,
p.bc=all_cond_BC$p, SIMPLIFY=FALSE)

2.2 Fonction EvalRepet

# Création de la fonction EvalRepet qui va appliquer les différents outils
# MIC sur les jeux de données générés juste au dessus
EvalRepet <- function(Re_Spectrum_data, group, groupN, nClust,
nLVPLSDA){
# Re_Spectrum_data : spectre prétraité
# group et groupN: classe et classe numérique des observations
# nClust : nombre de clusters pour le clustering
# nLVPLSDA: nombre de variables latentes pour le PLS-DA
# Inertia
Inertia.res = MBXUCL::Inertia(x = Re_Spectrum_data, y = groupN, print = FALSE)
# PCA
res.pca <- SVDforPCA(Re_Spectrum_data)
scores12 <- DrawScores(res.pca, type.obj = "PCA", color = group_num,
axes = c(1:2), main = "PCA scores")
scores34 <-DrawScores(res.pca, type.obj = "PCA", color = group, axes = c(3:4))
# clustering
ClustMIC.res = MBXUCL::ClustMIC(Intensities = Re_Spectrum_data, nClust = nClust,
Trcl = groupN, Dendr = FALSE)
# PLS-DA
PLSDA.res = PLSDA(x = Re_Spectrum_data, y = groupN,
nLV = nLVPLSDA, drawRMSEP = FALSE)
scoresPLS12 <- DrawScores(PLSDA.res, drawNames=TRUE, type.obj = "PLSDA",
createWindow=FALSE, main = paste0("PLSDA score plot for ", dataname),
color = group, pch = group,axes =c(1,2))
loadingsPLS12 <- DrawLoadings(PLSDA.res, type.obj = "PLSDA",
createWindow=FALSE, main = paste0("PLSDA loadings plot for", dataname),
axes = c(1:2), loadingstype="l", num.stacked = 2)
return(list(scores12=scores12,scores34=scores34,
Inertia.res=Inertia.res,
ClustMIC.res=ClustMIC.res,
PLSDA.res=PLSDA.res,
scoresPLS12=scoresPLS12,
loadingsPLS12=loadingsPLS12))
}
# Application de la fonction EvalRepet à toutes les matrices spectrales
# contenues dans res_Preprocessing et prétraitées avec des paramètres
# de correction de la baseline différents ---------------------
res.EvalRepet <- lapply(res_Preprocessing, FUN=EvalRepet, group = group,
groupN = group_num, nClust=4, nLVPLSDA=4)
names(res.EvalRepet) <- BCconditions
# Arrangement des sorties de la fonction EvalRepet ---------------------
# inertia
res.inertia <- sapply(res.EvalRepet, function(x) x[["Inertia.res"]][["Between_within"]]["Percentage",])
colnames(res.inertia) <- BCconditions
rownames(res.inertia) <- c("Inertia Between groups","Inertia Within groups", "Total inertia")
res.inertia <- round(res.inertia,2)
# PCA
PCAscores12 <- lapply(res.EvalRepet, function(x) x[["scores12"]])
names(PCAscores12) <- BCconditions
# clustering
res.clustering <- sapply(res.EvalRepet, function(x) x[["ClustMIC.res"]])
colnames(res.clustering) <- BCconditions
dimnam <- dimnames(res.clustering)
res.clustering <- matrix(data = unlist(res.clustering),
ncol = length(BCconditions), dimnames = dimnam,
byrow = FALSE)
res.clustering <- res.clustering[c(1,3,5,7),]
rownames(res.clustering) <- c("Dunn index", "Davies-Bouldin index", "Rand index", "Adjusted Rand index")
res.clustering <- round(res.clustering,2)
# PLS-DA
Q2_PLS <- sapply(res.EvalRepet, function(x) x[["PLSDA.res"]][["Q2"]])
Q2cum_PLS <- round(colMeans(Q2_PLS),2)
PLSscores12 <- lapply(res.EvalRepet, function(x) x[["scoresPLS12"]])
names(PLSscores12) <- BCconditions

2.3 Paramètres de la correction de la baseline \(\lambda\) et p

pander("lambda.bc (L) and p.bc (p) values for Baseline Correction")

lambda.bc (L) and p.bc (p) values for Baseline Correction

pander(all_cond_BC)
lambda p
1e+05 0.1
1e+06 0.1
1e+07 0.1
1e+08 0.1
1e+05 0.05
1e+06 0.05
1e+07 0.05
1e+08 0.05
1e+05 0.01
1e+06 0.01
1e+07 0.01
1e+08 0.01

2.4 Inertia

pander(res.inertia)
Table continues below
  L: 1e+05 p: 0.1 L: 1e+06 p: 0.1
Inertia Between groups 84.76 86.22
Inertia Within groups 15.24 13.78
Total inertia 100 100
Table continues below
  L: 1e+07 p: 0.1 L: 1e+08 p: 0.1
Inertia Between groups 87.49 89.2
Inertia Within groups 12.51 10.8
Total inertia 100 100
Table continues below
  L: 1e+05 p: 0.05 L: 1e+06 p: 0.05
Inertia Between groups 85.16 86.72
Inertia Within groups 14.84 13.28
Total inertia 100 100
Table continues below
  L: 1e+07 p: 0.05 L: 1e+08 p: 0.05
Inertia Between groups 88.26 89.82
Inertia Within groups 11.74 10.18
Total inertia 100 100
Table continues below
  L: 1e+05 p: 0.01 L: 1e+06 p: 0.01
Inertia Between groups 86.7 88.03
Inertia Within groups 13.3 11.97
Total inertia 100 100
  L: 1e+07 p: 0.01 L: 1e+08 p: 0.01
Inertia Between groups 89.54 90.07
Inertia Within groups 10.46 9.93
Total inertia 100 100
par(mar=c(2,6,2,2))
barplot(res.inertia["Inertia Between groups",],
main= "Inertia Between groups", col="blue", ylim = c(60,100), las=2)

2.5 PCA scores plots

pander("PCA scores plots")

PCA scores plots

PCAscores12
## $`L: 1e+05 p: 0.1`

## 
## $`L: 1e+06 p: 0.1`

## 
## $`L: 1e+07 p: 0.1`

## 
## $`L: 1e+08 p: 0.1`

## 
## $`L: 1e+05 p: 0.05`

## 
## $`L: 1e+06 p: 0.05`

## 
## $`L: 1e+07 p: 0.05`

## 
## $`L: 1e+08 p: 0.05`

## 
## $`L: 1e+05 p: 0.01`

## 
## $`L: 1e+06 p: 0.01`

## 
## $`L: 1e+07 p: 0.01`

## 
## $`L: 1e+08 p: 0.01`

2.6 Clustering

pander("clustering")

clustering

pander(res.clustering)
Table continues below
  L: 1e+05 p: 0.1 L: 1e+06 p: 0.1 L: 1e+07 p: 0.1
Dunn index 0.5 0.59 0.64
Davies-Bouldin index 1.01 0.92 0.8
Rand index 0.97 1 1
Adjusted Rand index 0.91 1 1
Table continues below
  L: 1e+08 p: 0.1 L: 1e+05 p: 0.05
Dunn index 0.73 0.49
Davies-Bouldin index 0.64 0.96
Rand index 1 0.97
Adjusted Rand index 1 0.91
Table continues below
  L: 1e+06 p: 0.05 L: 1e+07 p: 0.05
Dunn index 0.64 0.68
Davies-Bouldin index 0.84 0.76
Rand index 1 1
Adjusted Rand index 1 1
Table continues below
  L: 1e+08 p: 0.05 L: 1e+05 p: 0.01
Dunn index 0.89 0.68
Davies-Bouldin index 0.7 0.85
Rand index 1 1
Adjusted Rand index 1 1
Table continues below
  L: 1e+06 p: 0.01 L: 1e+07 p: 0.01
Dunn index 0.68 0.77
Davies-Bouldin index 0.69 0.7
Rand index 1 1
Adjusted Rand index 1 1
  L: 1e+08 p: 0.01
Dunn index 0.83
Davies-Bouldin index 0.67
Rand index 1
Adjusted Rand index 1
par(mar=c(6,2,2,2))
matplot(1:length(BCconditions),t(res.clustering), ylim=c(0,1),type="b",col=c("blue", "red", "forestgreen", "deepskyblue"),
pch=1:4, xaxt = "n", main = "Indices de clustering", xlab="")
legend("bottomleft", inset=0.01, legend=rownames(res.clustering),
pch=1:4, col=c("blue", "red","forestgreen", "deepskyblue"))
axis(side = 1, at = 1:length(BCconditions), labels = BCconditions,
las=2,cex.axis=0.8)

2.7 PLS-DA

pander("PLS-DA Q2cum")

PLS-DA Q2cum

pander(Q2cum_PLS)
Table continues below
L: 1e+05 p: 0.1 L: 1e+06 p: 0.1 L: 1e+07 p: 0.1 L: 1e+08 p: 0.1
0.93 0.94 0.93 0.95
Table continues below
L: 1e+05 p: 0.05 L: 1e+06 p: 0.05 L: 1e+07 p: 0.05 L: 1e+08 p: 0.05
0.94 0.93 0.95 0.94
L: 1e+05 p: 0.01 L: 1e+06 p: 0.01 L: 1e+07 p: 0.01 L: 1e+08 p: 0.01
0.95 0.94 0.94 0.94
par(mar=c(6,2,2,2))
barplot(Q2cum_PLS, main = "Q2 moyen pour toutes les réponses", col="blue", las=2)

pander("PLS-DA scores plots")

PLS-DA scores plots

PLSscores12
## $`L: 1e+05 p: 0.1`

## 
## $`L: 1e+06 p: 0.1`

## 
## $`L: 1e+07 p: 0.1`

## 
## $`L: 1e+08 p: 0.1`

## 
## $`L: 1e+05 p: 0.05`

## 
## $`L: 1e+06 p: 0.05`

## 
## $`L: 1e+07 p: 0.05`

## 
## $`L: 1e+08 p: 0.05`

## 
## $`L: 1e+05 p: 0.01`

## 
## $`L: 1e+06 p: 0.01`

## 
## $`L: 1e+07 p: 0.01`

## 
## $`L: 1e+08 p: 0.01`